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LONG-TERM  GOAL 

The  long-term  goal  of  this  project  is  to  develop  software  to  invert  acoustic  data  from  a  towed 
source/receiver  array  for  simultaneous  velocity  analysis,  detection  and  material  characterization  of 
small-scale  (7  to  15  cm)  scatterers  in  the  shallow  ocean  and  seabed  sediments.  Detection  will  be 
achieved  by  providing  a  reflector  map  of  the  target  region.  Material  characterization  is  provided 
through  estimates  of  angularly  dependent  reflection  coefficients  at  the  surface  of  the  target  over 
multiple  angles;  angle  versus  offset  analysis  provides  a  basis  for  characterizing  the  internal  material 
parameters  of  the  target. 

OBJECTIVES 

The  immediate  objective  is  to  complete  the  development  of  a  computer  code  for  towed  array  data  that 
back  propagates  the  data  with  a  wavespeed  that  varies  in  all  three  spatial  variables.  This  software  is 
also  to  be  integrated  with  software  to  carry  out  velocity  analysis  so  that  the  change  in  wavespeed  from 
the  water  column  to  the  seabed  can  be  properly  treated.  The  more  long-term  objectives  are  to  continue 
to  incorporate  developments  in  the  underlying  theory  into  the  software  being  developed  for  this 
application.  For  example,  the  theory  is  currently  being  extended  to  multi-pathing — the  case  in  which 
the  ray  trajectories  from  a  point  scatterer  at  depth  form  caustics  in  the  subsurface,  leading  to  multiple 
arrivals  at  the  same  observation  point  on  the  receiver  array. 

APPROACH 

The  particular  experiment  that  is  proposed  requires  the  acquisition  of  acoustic  responses  along  a  towed 
array  (about  10m)  of  sonophones  from  a  single  point  source.  The  experiment  is  to  be  repeated  on  a 
regular  near-surface  areal  grid.  Our  methods  are  based  on  a  high  frequency  asymptotic  inversion 
technique  that  is  called  Kirchhoff  inversion.  This  method  back  projects  the  observed  data  through 
integration  (summation)  over  diffraction  traveltime  curves  that  pick  out  the  observed  data  on  different 
traces  at  the  indicated  traveltimes.  The  distinctive  feature  of  our  approach  is  that  we  obtain  an  estimate 
of  the  specular  reflection  coefficient  and  specular  incidence  angle  from  the  peak  amplitude(s)  of  the 
output.  These  estimates  are  model-consistent  (back-projection  wavespeed)  and  require  amplitude 
control  of  the  source  signature.  On  other  scales,  this  method  has  achieved  great  success  in  seismic 
exploration,  with  lengths  two-to-three  orders  of  magnitude  higher-  and  frequencies,  correspondingly, 
two-to-three  orders  of  magnitude  lower-  than  those  used  for  the  experiments  here.  The  method  also 
has  application  to  flaw  characterization  in  solids — nondestructive  testing  and  evaluation. 
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WORK  COMPLETED 


We  have  previously  delivered  software  to  process  such  data  under  the  assumption  that  the  shallow 
water  and  seabed  environment  was  a  medium  with  a  wavespeed  that  was  only  depth  dependence;  first  a 
2.5D  code  was  developed,  then  a  3D  code.  The  designation,  2.5D,  refers  to  three  dimensional  wave 
propagation  over  an  ocean  and  seabed  model  that  is  assumed  to  have  two-dimensional  variability,  only. 
This  technique  allows  the  processing  of  a  single  line  of  data,  producing  a  hybrid  output  that  still  honors 
three-dimensional  propagation  properties— in  particular,  three-dimensional  geometrical  spreading  of  the 
amplitude.  The  current  software  uses  a  full  areal  survey  and  allows  for  three-dimensional  variation  in 
the  background  wavespeed.  It  also  will  produce  three-dimensional  images  of  the  target  scatterers  and 
the  seabed.  For  this  case,  the  first  challenge  was  to  develop  forward  modeling  software  that  produced 
sufficiently  accurate  amplitudes  as  well  as  traveltimes  in  the  forward  propagating  direction.  An 
additional  constraint  for  implementation  is  that  the  CPU  time  for  forward  modeling  should  be 
acceptable.  The  current  code  is  highly  accurate  in  both  traveltime  and  phase  for  forward  propagation, 
what  is  needed  for  the  inverse  problem.  By  smoothing  the  entire  model,  including  known  interfaces,  a 
code  running  in  acceptable  CPU  times  was  developed.  We  are  currently  experimenting  with  two 
methods  for  re-incorporate  transmission  coefficients  into  the  modeling  algorithm. 

A  crucial  ingredient  in  keeping  the  CPU  under  control  was  the  use  of  the  so-called  wavefront 
construction  method  for  modeling  the  traveltimes.  In  this  method,  an  appropriately  dense  set  of  rays  is 
set  off  near  the  source.  Triangular  tiles  on  the  wavefront  are  constructed  from  nearby  rays.  As  the 
wave  propagates  forward,  the  sizes  of  the  tiles  and  the  opening  angles  between  the  edges  of  the  tiles  are 
recorded.  When  these  numbers  reach  prescribed  thresholds,  new  rays  are  added.  In  this  manner,  rays 
are  only  added  as  they  are  needed  and  the  CPU  for  an  adequate  output  model  is  kept  to  a  minimum. 

RESULTS 

We  show  here  an  example  of  this  forward  modeling  code  applied  to  a  section  of  the  SEG/EAGE 
saltdome  model,  which  is  based  on  a  typical  US  Gulf  Coast  salt  structure.  The  complexity  of  the  3D 
salt  structure  makes  it  a  benchmark  model  used  to  evaluate  various  3D  forward  modeling  and  imaging 
algorithms.  Scaling  the  distances  down  by  a  factor  of  100  and  scaling  the  frequencies  up  by  the  same 
factor  provides  an  example  of  a  hard  structure  in  the  near  seabed.  Figure  1  shows  a  set  of  rays  from  a 
point  source.  The  velocity  in  the  salt  is  about  4482  m/s,  shown  in  black.  The  velocity  of  the  host 
medium  increases  with  depth,  but  is  much  lower  than  the  velocity  in  the  salt.  The  front  face  of  this 
figure  shows  a  vertical  slice  at  x  =  3.5km,  which  is  1.5km  from  the  source.  It  passes  through  the  salt 
structure.  The  rays  were  shot  off  with  a  uniform  distribution  of  take-off  angles.  The  feature  of  adding 
rays  to  maintain  ray  density  was  suppressed  in  this  figure,  so  that  we  could  see  the  relative  distribution 
of  the  set  of  rays  around  the  salt  structure.  Notice  that  the  density  of  the  post-critical  rays  that  do  not 
enter  the  salt  is  much  higher  than  that  of  the  rays  that  enter  the  salt  and  emerge  through  the  front  face. 
This  indicates  that  we  should  expect  the  amplitude  in  the  salt  to  be  lower  than  the  amplitude  right 
above  the  salt.  The  waves  that  propagate  through  the  salt  generate  first  arrivals  at  near-salt  and  sub-salt 
regions.  These  waves  dominate  the  solution  of  first  arrival  eikonal  solvers. 

Figure  2  is  one  wavefront  recorded  at  0.9  s  with  the  source  located  at  the  center  of  the  upper  surface. 

Its  lower  part  changes  from  the  typical  spherical  shape,  due  to  the  presence  of  the  high  velocity  salt. 
This  structure  causes  portions  of  the  wavefront  to  bulge  and  expand,  especially  those  parts  that  have 


penetrated  the  salt.  The  boundaries  of  this  bulge  are  dominated  by  the  pseudo-head-wave  arrivals, 
which  continuously  connect  the  portions  of  the  wavefront  that  penetrated  the  salt  with  those  that 
traveled  only  through  the  sediments.  However,  head-waves  are  low-energy  arrivals  that  are  not  useful 
for  imaging  applications.  Waves  traveling  directly  from  the  source,  without  going  through  the  salt, 
carry  more  energy  than  rays  that  propagate  through  the  salt. 

Figure  3  is  an  amplitude  plot  on  the  same  vertical  slice  as  the  front  face  of  Figure  1.  Here,  the  red  and 
yellow  colors  indicate  high  amplitude  and  the  purple  color,  low  amplitude.  The  high  and  low 
amplitude  zones  agree  with  what  we  expected  from  the  ray  plot  in  Figure  1.  From  examples  provided 
in  earlier  reports  to  ONR,  we  know  that  the  amplitudes  near  turning  are  not  as  accurate  as  the 
transmitted  amplitudes,  but  they  are  close  enough  qualitatively  to  confirm  the  agreement  we  see 
between  these  two  figures. 

As  noted  above,  even  the  reflectors  are  smoothed  in  preparation  for  the  ray  tracing  in  order  to  generate 
stable  amplitudes.  Thus,  it  is  necessary  to  insert  transmission  and  reflection  coefficients  in  a  non- 
traditional  manner.  Here,  we  use  a  scheme  that  updates  the  transmission  coefficients  pointwise.  That 
is,  an  approximate  transmission  coefficient  is  introduced  at  every  step  of  the  Runge-Kutta  routine  used 
for  the  ray  tracing.  Figure  4  depicts  level  curves  of  amplitude  for  a  simple  horizontal  interface  model, 
velocity,  2000m/s  above,  2500m/s,  below.  The  amplitude  is  generated  from  ray  theory  for  the 
discontinuous  model  with  level  curves  shown  as  solid  curves.  Then,  the  numerical  output  for  the 
smoothed  velocity  model  is  generated,  with  level  curves  depicted  with  dashed  curves.  The  numbers 
are  the  negative-log-amplitude.  In  this  case,  the  transmission  coefficient  is  equal  to  one.  Thus,  the 
method  compares  the  close  fit  and  stability  of  the  ray-theory  Jacobian  calculations  for  the  smoothed 
and  discontinuous  models.  They  are  seen  to  agree  extremely  well  in  the  forward  propagation  direction- 
— exactly  the  direction  of  interest  in  the  ultimate  application  to  inversion.  Figure  5  depicts  the 
amplitudes  again.  However,  now,  the  transmission  coefficient  is  included  in  the  analytically  computed 
result  and  the  pointwise  generation  of  the  transmission  coefficient  is  included  in  the  numerical 
modeling  result  produced  by  our  code.  Again,  the  agreement  in  the  forward  modeling  direction  is 
extremely  good.  Finally,  we  show  a  traveltime  comparison  in  Figure  6.  We  can  see  the  head  wave 
produced  by  the  analytic  solution  (solid  curve)  and  the  approximate  head  wave  produced  by  the 
computer  code  (dashed  curve).  Except  for  the  post-critical  region,  all  of  the  amplitude  and  traveltime 
comparisons  are  accurate  enough  for  proposed  use  in  the  inversion  code. 

IMPACT/APPLICATION 

Inversion  of  towed  array  data  under  the  assumption  of  a  fully  heterogeneous  background  wavespeed 
has  resided  in  commercial  processing  companies  servicing  the  oil  industry  and  in  the  oil  industry, 
itself.  We  know  from  our  connections  to  the  oil  industry  that  the  forward  modeling  program  that  we 
are  developing  is  competitive  with  modeling  programs  available  in  the  industry.  Such  forward 
modeling  programs  are  the  linchpins  on  which  one  builds  a  corresponding  inversion  code.  Thus,  we 
are  developing  code  that  will  be  generally  available  as  free  software  through  the  Center  for  Wave 
Phenomena  web  site.  We  are  aware  of  industrial  codes  available  in  the  industry  that  have  the  imaging 
capability  of  the  code  that  we  are  developing,  but  do  not  have  the  parameter  estimation  capability  of 
our  code. 


TRANSITIONS 


Our  ties  to  the  oil  industry  through  a  consortium  of  30  companies  provide  a  direct  line  for  technology 
transfer  and  prompt  and  extensive  testing  of  our  science  and  its  computer  implementation.  We  know 
that  our  methods  have  led  to  industry-standard  codes  for  processing  seismic  exploration  data  to  aid  in 
the  identification  of  fluid  traps  in  the  Earth.  Confirmation  of  the  efficacy  of  our  methods  comes 
through  annual  renewal  by  our  sponsors. 

RELATED  PROJECTS 

In  this  program,  we  assume  that  the  background  wavespeed  is  known.  Under  other  ONR  support,  we 
have  developed  a  three  dimensional  migration  velocity  analysis  method  to  help  determine  that 
background  wavespeed.  Also,  there  are  standard  preprocessing  techniques  designed  to  reduce  the  size 
and  complexity  of  the  data  as  well  as  to  carry  out  preliminary  analysis  in  advance  of  full  3D  data 
processing.  We  have  found  that  many  of  those  methods  fall  in  a  single  category  that  we  call  data 
mapping.  This  is  a  technique  for  transforming  data  from  an  input  source/receiver  configuration  and 
presumed  background  earth  model  to  a  different  source/receiver  configuration  and,  possibly,  a  different 
earth  model.  This  research  was  initiated  under  the  now-defunct  ACTEDOE  program.  It  is  continuing 
under  oil  industry  support  through  our  consortium  project.  We  are  developing  a  theory  that  we  call 
Kirchhoff  data  mapping  (KDM).  We  have  a  general  platform  for  this  process  and  have  developed 
implementations  in  the  following  list,  with  applicability  and  limitations  as  further  implied  by  the  titles 
in  the  references,  below. 

1 .  Mapping  of  data  gathered  at  finite,  fixed  offset  (common  offset,  bistatic)  between  source  and 
receiver  to  zero  offset  (monostatic)  data. 

2.  Downward  continuation  of  receivers  (or  sources). 

3.  Mapping  of  common  offset  data  to  common  shot  data.  This  output,  here  provides  synthetic 
data  as  if  the  receiver  array  were  as  long  as  the  entire  survey,  rather  than  as  long  as  the  given 
towed  array. 

Other  implementations  are  under  development  at  this  time.  There  is  high  interest  in  this  subject  in  the 
oil  industry. 
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Figure  1.  Rays  from  a  point  source 
propagating  through  a  salt  structure. 


Figure  4.  Amplitude  comparison  for  horizontal 
interface  model,  velocity,  2000m/s  (above), 
2500m/s  (below).  Transmission  coefficient  =  1. 
Solid  curve:  ray  theoretic  amplitude  for 
discontinuous  medium;  dashed  curve:  smoothed 
medium.  Numbers:  -log  amplitude 
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Figure  2.  A  wavefront  or  equal  travel  surface  Figure  5.  Amplitude  comparison  with  true  trans- 
at .  9  seconds.  mission  coefficient. 


Y  (km) 


X  (km) 

12  3  4 


Figure  3.  Amplitude  plot  on  the  same  front  face 
as  in  Figure  1. 


Figure  6.  Traveltime  comparison. 


